knitr::opts_chunk$set(echo = TRUE,
warning= FALSE,
message = FALSE,
error=FALSE)
library(knitr)
library(dplyr,quietly = T,warn.conflicts = F,verbose = F)
library(data.table,warn.conflicts = F, verbose=F)Warning: package 'data.table' was built under R version 4.3.3library(sf,quietly = T, verbose=F)Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUENous partons des fichiers détail de l’INSEE car, à la différence des tableaux prédéfinis, ils permettent virtuellement toutes les formes de croisement d’indicateurs. Ils sont évidemment très volumineux mais ici vous avez la chance de bénéficier d’une extraction des communes sur lesquelles vous allez travailler
Comme les fichiers sont d’assez petite taille, on se contentera d’utiliser la fonctionread.csv qui fait partie du R-base :
t1<-Sys.time()
tab<-read.csv("data/rplogmt_13055.csv")
t2<-Sys.time()
paste ("chargement effectué en",t2-t1,"secondes")
dim(tab)
class(tab)Dans l’exemple choisi (Marseille) le tableau ne comporte que 182735 individus et se charge en moins d’une seconde
Toutefois, si vous êtes amenés à traiter directement des fichiers de l’INSEE comportant plusieurs millions d’enregistrement vous aurez intérêt à utiliser des fonctions plus puissantes comme la fonction read_csv2 du package readr ou, miuex encore, la fonction fread du package data.table
Nos tableaux peuvent être enregistés au format interne de .R ce qui réduira considérablement leur taille par rapport au fichier texte au format .csv qui pèse en général beaucoup plus.
Comme ces données vont nous servir durant tout le projet, elles seront stockées dans le dossier res qui sera spécifique à l’exemple retenu.
## Sauvegarde au format RDS
saveRDS(object = tab,
file = "res/tabmen.RDS")Il ne faut surtout pas oublier le fichier des métadonnées qui va permettre de recoder facilement tous les facteurs et de décoder les chiffres correspondant aux classes. On va donc le transformer au format R puis l’enregistrer également dans le dossier res. Noter qu’on utilise cette fois-ci la fonction read.csv2car le séparateur de colonne n’est plus la virgule (ce qu’on lit avec read.csv) mais le point virgule.
# Lecture du fichier de métadonnées
meta<-read.csv2("meta/Dictionnaire variables_logemt_2020.csv")
# Enregistrement dans le dossier data
saveRDS(object = meta,
file = "res/tabmen_meta.RDS")Les contours des unités spatiales correspondant aux codes de l’INSEE sont produits par l’IGN et disponibles sur le site géoservice en accès libre. En général on récupère un très gros fichier pour la France entière ce qui demande un peu de travail pour extraire juste la zone d’étude. Ici vous avez de la chance puisqu’on vous a déjà extrait le fichier relatif à votre commune
La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets ubniques rassemblant à la fois
Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou dans d’autres formats standards comme GeoJson, la première tâche consiste donc à les convertir au formt sf afin de pouvoir les utiliser facilement dans R. L’importation se fait à l’aide de l’instruction st_read en indiquant juste le nom du fichier .shp à charger. Les autres fichiers (.dbf ou .proj) seront lus également et intégrés dans l’objet qui hérite de la double classe data.frame et sf.
library(sf)
map_iris <- st_read("data/iris_13055.gml")Reading layer `iris_13055' from data source
`/Users/claudegrasland1/worldregio/hackaton2024/data/iris_13055.gml'
using driver `GML'
Simple feature collection with 393 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 881250.9 ymin: 6232998 xmax: 905587.3 ymax: 6257715
Projected CRS: RGF93 Lambert 93dim(map_iris)[1] 393 8class(map_iris)[1] "sf" "data.frame"head(map_iris,2)Simple feature collection with 2 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 892700.6 ymin: 6245602 xmax: 893856.1 ymax: 6246686
Projected CRS: RGF93 Lambert 93
gml_id INSEE_COM NOM_COM IRIS CODE_IRIS
1 iris_13055.0 13206 Marseille 6e Arrondissement 301 132060301
2 iris_13055.1 13206 Marseille 6e Arrondissement 602 132060602
NOM_IRIS TYP_IRIS geometryProperty
1 Bergers H MULTIPOLYGON (((893856.1 62...
2 Vauban-Fénelon H MULTIPOLYGON (((892867.6 62...Pour faciliter la compatibilité avec d’autres applications, la 8e colonne qui contient le contour des zones est renommée “geometry” plutôt que “geometryProperty”.
st_geometry(map_iris) <-"geometry"Le nouveau tableau ne comporte plus que 393 unités spatiales et 8 colonnes On peut visualiser rapidement le résultat à l’aide de la fonction plotdu package sf :
par(mar=c(0,0,0,0))
plot(map_iris$geometry)
On sauvegarde le résultat dans notre dossier res au format interne de R :
saveRDS(object = map_iris,
file = "res/map_iris.RDS")Comme nous serons amenés à travailler à plusieurs échelles, nous produisons tout de suite un fonds de carte des arrondissements en utilisant les fonctions d’agrégation du packages sf combinées avec celles de dplyr.
map_arr <- map_iris %>% group_by(INSEE_COM) %>%
summarise(NOM_COM = min(NOM_COM)) %>%
st_as_sf()on vérifie que l’agrégation s’est bien passée :
par(mar=c(0,0,0,0))
plot(map_arr$geometry)
Et on sauvegarde le fonds de carte
saveRDS(object = map_arr,
file = "res/map_arr.RDS")Enfin, on construit un fonds de carte de l’ensemble de la commune étudiée
# Agrégation
map_com <- map_iris %>% mutate(tot = "Marseille") %>%
group_by(tot) %>%
summarise() %>%
st_as_sf()
# Visualisation
par(mar=c(0,0,0,0))
plot(map_com$geometry)
# Sauvegarde
saveRDS(object = map_com,
file = "res/map_com.RDS")La fonction plot du package sf permet de superposer des couches cartographiques avec l’instruction add=TRUE. Elle permet également de régler les couleurs des polygones (col=) ainsi que la couleur et l’épaisseur des frontières (border= et lwd=). Lorsqu’on ne souhaite pas remplir les polygones on indique col=NULL.
par(mar=c(0,0,0,0))
plot(map_iris$geometry, col= "lightyellow", border = "gray80", lwd=0.5)
plot(map_arr$geometry, col= NA, border = "gray50", lwd=1, add=TRUE)
plot(map_com$geometry, col= NA, border = "gray30", lwd=2, add=TRUE)
Les données sur le logement social sont été téléchargées sur le site du ministère du développement durable :
On charge le fichier qui comporte 182735 observation et 70 variables et on le sauvegarde
rpls<-read.csv("data/rpls_13055.csv")
head(rpls) X.1 DROIT DEPCOM CODEPOSTAL LIBCOM NUMVOIE INDREP
1 1 1 13055 13009 Marseille 9e Arrondissement <NA> <NA>
2 2 1 13055 13009 Marseille 9e Arrondissement <NA> <NA>
3 3 1 13055 13009 Marseille 9e Arrondissement <NA> <NA>
4 4 1 13055 13009 Marseille 9e Arrondissement <NA> <NA>
5 5 1 13055 13009 Marseille 9e Arrondissement <NA> <NA>
6 6 1 13055 13009 Marseille 9e Arrondissement <NA> <NA>
TYPVOIE NOMVOIE NUMAPPT NUMBOITE ESC COULOIR ETAGE
1 <NA> 55, Traverse de la Gouffonne 103 NA <NA> <NA> 0
2 <NA> 55, Traverse de la Gouffonne 113 NA <NA> <NA> 1
3 <NA> 55, Traverse de la Gouffonne 114 NA <NA> <NA> 1
4 <NA> 55, Traverse de la Gouffonne 202 NA <NA> <NA> 0
5 <NA> 55, Traverse de la Gouffonne 123 NA <NA> <NA> 2
6 <NA> 55, Traverse de la Gouffonne 124 NA <NA> <NA> 2
COMPLIDENT ENTREE BAT IMMEU COMPLGEO LIEUDIT QPV
1 <NA> <NA> <NA> 184-BT1-1 RESIDENCE LE MASSALIA I <NA> 2
2 <NA> <NA> <NA> 184-BT1-1 RESIDENCE LE MASSALIA I <NA> 2
3 <NA> <NA> <NA> 184-BT1-1 RESIDENCE LE MASSALIA I <NA> 2
4 <NA> <NA> <NA> 184-BT1-2 RESIDENCE LE MASSALIA I <NA> 2
5 <NA> <NA> <NA> 184-BT1-1 RESIDENCE LE MASSALIA I <NA> 2
6 <NA> <NA> <NA> 184-BT1-1 RESIDENCE LE MASSALIA I <NA> 2
TYPECONST NBPIECE SURFHAB CONSTRUCT LOCAT PATRIMOINE ORIGINE FINAN FINANAUTRE
1 C 3 65 1995 1998 1998 1 12 12
2 C 3 65 1995 1998 1998 1 12 12
3 C 5 100 1995 1998 1998 1 12 12
4 C 3 65 1995 1998 1998 1 12 12
5 C 3 69 1995 1998 1998 1 12 12
6 C 3 69 1995 1998 1998 1 12 12
CONV NUMCONV DATCONV NEWLOGT CUS DPEDATE DPEENERGIE
1 1 13204199785123110750722399 NA NA NA NA B
2 1 13204199785123110750722399 NA NA NA NA B
3 1 13204199785123110750722399 NA NA NA NA B
4 1 13204199785123110750722399 NA NA NA NA B
5 1 13204199785123110750722399 NA NA NA NA B
6 1 13204199785123110750722399 NA NA NA NA B
DPESERRE SRU_EXPIR SRU_ALINEA CODSEGPATRIM LIBSEGPATRIM PMR REG
1 C 2031 NA SEG001 Marseille Prov. Metropol 1 93
2 C 2031 NA SEG001 Marseille Prov. Metropol 1 93
3 C 2031 NA SEG001 Marseille Prov. Metropol 1 93
4 C 2031 NA SEG001 Marseille Prov. Metropol 1 93
5 C 2031 NA SEG001 Marseille Prov. Metropol 1 93
6 C 2031 NA SEG001 Marseille Prov. Metropol 1 93
LIBREG DEP LIBDEP EPCI
1 Provence-Alpes-Côte d'Azur 13 Bouches-du-Rhône 200054807
2 Provence-Alpes-Côte d'Azur 13 Bouches-du-Rhône 200054807
3 Provence-Alpes-Côte d'Azur 13 Bouches-du-Rhône 200054807
4 Provence-Alpes-Côte d'Azur 13 Bouches-du-Rhône 200054807
5 Provence-Alpes-Côte d'Azur 13 Bouches-du-Rhône 200054807
6 Provence-Alpes-Côte d'Azur 13 Bouches-du-Rhône 200054807
LIBEPCI EPSG X Y PLG_QP PLG_IRIS2021
1 Métropole d'Aix-Marseille-Provence 2154 895698.9 6241946 HZ 404
2 Métropole d'Aix-Marseille-Provence 2154 895698.9 6241946 HZ 404
3 Métropole d'Aix-Marseille-Provence 2154 895698.9 6241946 HZ 404
4 Métropole d'Aix-Marseille-Provence 2154 895698.9 6241946 HZ 404
5 Métropole d'Aix-Marseille-Provence 2154 895698.9 6241946 HZ 404
6 Métropole d'Aix-Marseille-Provence 2154 895698.9 6241946 HZ 404
PLG_ZUS PLG_ZFU PLG_QVA PLG_CODE_COMMUNE2021 QUALITE_VOIE QUALITE_NUMERO
1 HZ CSZ HZ 13209 NA NA
2 HZ CSZ HZ 13209 NA NA
3 HZ CSZ HZ 13209 NA NA
4 HZ CSZ HZ 13209 NA NA
5 HZ CSZ HZ 13209 NA NA
6 HZ CSZ HZ 13209 NA NA
QUALITE_XY DISTANCE_PRECISION QUALITE_QP QUALITE_IRIS QUALITE_ZUS QUALITE_ZFU
1 11 <NA> 1 1 1 NA
2 11 <NA> 1 1 1 NA
3 11 <NA> 1 1 1 NA
4 11 <NA> 1 1 1 NA
5 11 <NA> 1 1 1 NA
6 11 <NA> 1 1 1 NA
QUALITE_QVA COMAQP COMRIL
1 1 1 1
2 1 1 1
3 1 1 1
4 1 1 1
5 1 1 1
6 1 1 1saveRDS(rpls,"res/rpls.RDS")Le fichier comporte deux colonnes X et Y dont on ignore la projection. On va tester leur superposition sur les fonds de carte précédents pour voir s’il s’agit de la même.
par(mar=c(0,0,0,0))
plot(map_iris$geometry, col= "lightyellow", border = "gray80", lwd=0.5)
plot(map_arr$geometry, col= NA, border = "gray50", lwd=1, add=TRUE)
plot(map_com$geometry, col= NA, border = "gray30", lwd=2, add=TRUE)
points(rpls$X, rpls$Y, col="red",pch=20, cex=0.2)
Bingo, la superposition est parfaite donc la projection est la même. Il s’agit de la Lambert93 dont le CRS peut-être affiché par la commande st_crs() du package sf :
st_crs(map_iris)Coordinate Reference System:
User input: RGF93 Lambert 93
wkt:
PROJCRS["RGF93 Lambert 93",
BASEGEOGCRS["RGF93 geographiques (dms)",
DATUM["Reseau Geodesique Francais 1993 v1",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["IGNF","RGF93G"]],
CONVERSION["LAMBERT-93",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",46.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",3,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",44,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",49,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",700000,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",6600000,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["NATIONALE A CARACTERE LEGAL"],
AREA["FRANCE METROPOLITAINE (CORSE COMPRISE)"],
BBOX[41,-5.5,52,10]],
ID["IGNF","LAMB93"]]On va alors pouvoir créer un fonds de carte géométrique de ces logements sociaux en indiquant la projection du fonds de carte des Iris.
# Transformation du data.frame en objet sf
rpls_geom<-st_as_sf(rpls,coords = c("X","Y"))
# Récupération du code EPSG de la projection
st_crs(rpls_geom)<-st_crs(map_iris)
# Sauvegarde au format .RDS
saveRDS(rpls_geom, "res/rpls_geom.RDS")Nous avons désormais un dossier res qui comporte tous les éléments utiles pour l’analyse de notre zone d’étude.
list.files("res")[1] "map_arr.RDS" "map_com.RDS" "map_iris_hlm.RDS" "map_iris.RDS"
[5] "rpls_geom.RDS" "rpls.RDS" "tabmen_meta.RDS" "tabmen.RDS"